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We present the theoretical study on non-equilibrium (NEQ) fluctuations for diffusion dynamics 
in high dimensions driven by a linear drift force. We consider a general situation in which NEQ is 
caused by two conditions: (i) drift force not derivable from a potential function and (ii) diffusion 
matrix not proportional to the unit matrix, implying non-identical and correlated multi-dimensional 
noise. The former is a well-known NEQ source and the latter can be realized in the presence of 
multiple heat reservoirs or multiple noise sources. We develop a statistical mechanical theory based 
on generalized thermodynamic quantities such as energy, work, and heat. The NEQ fluctuation 
theorems are reproduced successfully. We also find the time-dependent probability distribution 
function exactly as well as the NEQ work production distribution P{W) in terms of solutions of 
nonlinear differential equations. In addition, we compute low-order cumulants of the NEQ work 
production explicitly. In two dimensions, we carry out numerical simulations to check out our 
analytic results and also to get P(VV). We find an interesting dynamic phase transition in the 
exponential tail shape of P(VV), associated with a singularity found in solutions of the nonlinear 
differential equation. Finally, we discuss possible realizations in experiments. 

PACS numbers: 05.70. Ln, 05.10.Gg, 05.40.-a 



I. INTRODUCTION 

There have been great interests in non-equilibrium 
(NEQ) statistical mechanics for last decades since the 
discovery of the fluctuation theorem for entropy produc- 
tion. The first discovery was made on a deterministic 
NEQ dynamics governed by the SLLOD equation [ll-Q . 
Later on, the fluctuation theorems of various types were 
found to be a universal feature for a wide class of NEQ 
systems, which are governed by both deterministic 0, Hj 
and stochastic dynamics (gI-ITH . Jarzynski found an inter- 
esting relation between NEQ work and equilibrium free 
energy , which was later proved to be a special case of 
Crooks fluctuation theorems 0, [3] ■ Since then, a num- 
ber of stimulated studies have been published up to now 
regarding the fluctuation theorems and related phenom- 
ena [Il-i^. 

The diffusion dynamics is distinguished from the Brow- 
nian dynamics. The former has only position-like state 
variables that have even parity under the time reversal, 
while the latter has pairs of position and momentum with 
even and odd parities respectively. In this work, we con- 
sider the diffusion dynamics with two important condi- 
tions which drive the system into a NEQ steady state 
(NESS): (i) drift force not derivable from a potential 
function (ii) non-identical and correlated noise. These 
two conditions can be realized only in high dimensions 
(not possible in one dimension). Unusual results were 
reported for the dynamics with the combination of these 



* lckwon@mju. ac.kr] 
t jdnoh@uos.ac! kr| 
t |hgpark@kias.re.kr| 



two conditions |21[ - |23| . In particular it was found that 
the zero mass limit and the over-damping limit are dif- 
ferent in reducing the Kramers equation to the Fokker- 
Planck equation. 



Kwon, Ao and Thouless [2J| studied the diffusion dy- 
namics with a linear drift force in high dimensions and 
found the probability distribution function (PDF) for the 
NESS exactly. They found that a circulating probabil- 
ity current can exist at the steady state, violating the 
detailed balance. One example is a non-zero torque gen- 
erated in a nano heat engine in contact with two dif- 



ferent heat reservoirs [25[. For a nonlinear drift force, 



it has recently been found via a perturbation theory 
that there exists an additional current, that is absent 
for the linear case, due to the combination of the force 
non-linearity and the multi-dimensional noise correlation 
[26j . It moves the probability maximum away from the 
fixed point at which the force is zero. This novel current 
has the same origin with the noise- induced current, trans- 
porting drugs or molecules in biologicalsystems, studied 
by Prost et al. and Doering et al. . 

We revisit the diffusion dynamics with a linear drift 
force from the point of view of the NEQ fluctuation the- 
orem, as one of a few analytically solvable cases far from 
equilibrium. We consider the Langevin equation 



(1) 



where q = {qi,q2, ■ ■ ■ ,qd)^ is a state vector in d dimen- 
sions with the superscript T denoting the transpose of 
a given vector or matrix. We restrict ourselves to the 
case in which q has even parity under the time rever- 
sal, i.e., there are no momentum-like variables. ^ = 
(?iiC2j ■ ■ ■ i^d)"^ is a white noise vector with zero mean 



satisfying {^{t)^^ {t')) = 2D5{t - t') where 



is the 
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noise average. D is defined as a d x d diffusion matrix 
that is symmetric, positive definite, and g- independent. 
The first NEQ condition (i) is given by tlie drift force 
/ 7^ —V^{q) wfiere $ is a scalar function of q. The con- 
dition (ii) leads to the case where the diffusion matrix is 
not proportional to the unit matrix, D I, va contrast 
to the conventional thermal noise with D cx (ksT)!. 

This Langevin equation describes a general stochastic 
system far from equilibrium without conventional energy 
or temperature. So we first define generalized thermody- 
namic quantities; energy, work, and heat properly. With 
these definitions, we successfully reproduce the NEQ 
fluctuation theorems. We can also get analytic expres- 
sions for many interesting quantities such as the time- 
dependent PDF P{q,t), the two-time correlation func- 
tions, and the NEQ work production distribution P{W). 

The Jarzynski equality can be shown directly and the 
cumulants for the NEQ work production W are calcu- 
lated explicitly up to the second order. More interest- 
ingly, we find the exponential tail shape of P{W) with a 
power-law prefactor, which undergoes a dynamic phase 
transition as the time increases. This phase transition 
turns out to be associated with a singularity of solutions 
of a nonlinear differential equation (NLDE). We solve 
this NLDE numerically to reveal the details of the expo- 
nential tail shape of P{W) and also its dynamic phase 
transition. 

This paper is organized as follows. In Sec. II, the gen- 
eralized thermodynamic quantities are defined with the 
corresponding fluctuation theorems. In Sec. Ill, we ob- 
tain P{q, t) and the two-time correlation functions ex- 
actly. In Sec. IV, we derive the analytic expression for the 
generating function of P{W) in terms of solutions of the 
NLDE. In Sec. V, we calculate the cumulants of the work 
production. In Sec. VI, we take the two-dimensional dif- 
fusion dynamics as a simple example and calculate the 
generating function by solving the NLDE numerically. 
The dynamic phase transition of the tail shape of PiW) 
is discussed. We also present the results from the direct 
numerical integration of the Langevin equation, which 
agree with the analytical results. Finally, in Sec. VII, we 
summarize our results and discuss novel features found 
for the linear diffusion dynamics in high dimensions and 
the possibility of realization in experiments. 



II. FLUCTUATION THEOREMS FOR THE 
DIFFUSION DYNAMICS IN HIGH DIMENSIONS 

We consider the Fokker-Planck equation for the diffu- 
sion dynamics associated with the Langevin equation in 
Eq. (P, 

^^^=y-{~f{q) + D-V)P{q,t) , (2) 

where we use the dot notation for the product of a vec- 
tor and a matrix, or between vectors, but not between 



matrices. For example, we write q^ ■ AB ■ q where A, B 
are matrices. 

Writing the steady state solution as Pst oc e~*=*^''\ 
we define equilibrium as the steady state satisfying the 
detailed balance such that 

n(g',t';<z,t)e-*-(') =n(q,t;q',i')e-*-(«'', (3) 

where Ii{q\t']q,t) is the conditional probability for the 
transition from state q at time t to state q' at time t' . 

We can show that the necessary and sufficient condi- 
tion for the detailed balance reads as 

f(q) = -D-V^st{q). (4) 

i.e., the vanishing probability current J = { f—D-'S/)Pst = 
at the steady state. Defining the force matrix F as 

F^p - -V/j/„, (5) 

with = we can rewrite the detailed balance con- 
dition as 

FD = (FDf = DF^ , (6) 

where VaV^<I>st = V^Vc^st is used. 

The detailed balance condition is always satisfied in 
one dimension. In higher dimensions, however, this does 
not hold in general due to two possible sources. One is 
the asymmetry oi F {F ^ F'^)i which happens when / 
is not derivable from a scalar potential. Another comes 
from the diffusion matrix which is not proportional to the 
unit matrix, so the noises are not identical in components 
(Daa 7^ Dpi3 for a ^ j3) and also may be correlated 
(D^p + 0). 

We note that the detailed balance condition is satis- 
fied (subsequently, equilibrium can be achieved) not only 
with symmetric F and D (x I, but also with the specific 
combination of general F and D I satisfying Eq. ([6]). 
In equilibrium, the steady-state distribution should be 
given by the Boltzmann distribution, so can be inter- 
preted as energy (we set ksT = I for convenience). Then 
it is natural to define a generalized force as the negative 
derivative of the energy function, i.e. ■ f = — V<f>st 
from Eq. ([5]), especially for D I. 

In general, the detailed balance condition is not satis- 
fied and the system is driven into a NESS. In this case, 
Eq. Q should be modified as 

■ f ^ -V^q) + g{q) , (7) 

with nonzero g(q) which can not be derivable from a 
scalar function {g{q) ^ — V^((j')). Wc interpret g as 
a generalized NEQ driving force. However, there is no 
unique way of defining the NEQ driving force g as well 
as the energy-like function $. As the NESS needs not 
be governed by the Boltzmann distribution, the energy 
function $((7) does not have to be the same as $^4(9) 
in general. In fact, one may choose an arbitrary $((7) 
and g(g) accordingly in the system described only by the 
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stochastic equations like Eqs. ([T|) and We will come 
back to this issue later. 

Following the path integral formalism of Onsager and 
Machlup [2a| , the conditional probability that the system 
evolves along a path ^(t) for < r < t starting from an 
initial state q{0) is given by 

The time-reverse path is given by ^(t) = q(t — r) with 
the initial state g(0) = q{t). Then the ratio of the con- 
ditional probability for the forward path q(r) to that for 
the reverse path ^(t) is found as 

mr)Yo 

^ g-A4.+W[?] ^ g-Q[g] ^ ^g^ 

where A$ ~ ^{q{t)) — ^{q{0)) is the energy difference 
between at the final and initial time. We interpret Wlq] 
as the work production done by the NEQ force g along 
the path g(r), 

WM = £rfr(^g^..g+^^ . (10) 

The second term appears only when $ has an explicit 
time dependence, which is the case of the Jarzynski type 
work |12| . 

Accordingly we define Q = A$ — W and interpret it 
as the heat transferred into the system from the reservoir 
along the path. We can rewrite 

Q[q]= f dT{~f -D-^ -q + f -C) (11) 
Jo 

where C = -{q- ,f) ^ D-^ ■ This definition of heat 
can be understood in terms of the corresponding Brown- 
ian dynamics. The diffusion dynamics can be regarded as 
the Brownian dynamics in the overdamped limit (or zero 
inertia limit). In this case, the corresponding Brownian 
dynamics is given by 

mq^^j-q + D-'-f + C, (12) 

where the friction matrix ^ ~ D^^ and the diffusion ma- 
trix for the noise C = D~^£^ becomes = D~^DD^^ = 
D~^. Therefore, the generalized Einstein relation is sat- 
isfied (£)(; = ksTj with our presetting of /cbT = 1). As 
noticed, ■ f plays the role of force. From the point 
of view of the Brownian dynamics, C — • ^ is a ran- 
dom diffusive force exerted by noise. Then Q is the sum 
of dissipative work due to the random collisions and the 
diffusive work due to noise, known as the Onsager heat. 
Therefore the interpretation of Q as the heat production 
is reasonable. 

Up to now, we introduced the appropriate definitions 
of generalized energy, work, and heat for an arbitrary 
diffusion dynamics without conventional thermodynamic 



quantities. Note that the work and the heat produc- 
tion depend on the path, i.e., functionals of the path, 
while the energy is a state function and independent of 
the path. One may restore the temperature by scaling 
the diffusion matrix as £) — (fc^T)!), then ^st scales as 
^st/{kBT). ksT can parametrize the overall strength of 
noise and can be interpreted as an effective temperature, 
if necessary. There are other temperature-like parame- 
ters, the eigenvalues of D, which may be different each 
other. The presence of multi and heterogeneous temper- 
atures might also cause NEQ, though not always. 

Eq. ^ is the fundamental equation from which vari- 
ous fluctuation theorems can be derived. It replaces the 
detailed balance relation for equilibrium, so is referred 
to as the detailed fluctuation relation or the generalized 
detailed balance relation for NEQ. If we choose an initial 
PDF cx e^**^'^ (Boltzmann distribution with the energy 
function <&((?)), we can easily show the Crooks fluctuation 
theorem @: 

{A[q])F = a[g]e^[^1)fl e-^^ , (13) 

where A[q] is any functional of the path q^r) and A[q\ = 
A[q]. The free energy is deflned as e~'^ = / dqe'"^ and 
Eq. Uni) guarantees that W[q] ^ -W[q]. The subscripts 
F and R denote the averages along the forward q{T) and 
the reverse g(r) path respectively. The free energy dif- 
ference A J" = J'(i) — J-'(O) does not vanish only when <I> 
has an explicit time dependence (^ 7^ 0). In this paper, 
we assume no explicit time dependence in / and I?, so 
we always find AT = 0. 

III. TRANSIENT STATE FAR FROM STEADY 
STATE FOR LINEAR DIFFUSION DYNAMICS 

From now on, we focus on the case of a linear drift 
force, 

f = -F-q, (14) 

which can be analytically tractable. The force matrix F 
is constant of state q and time t. The exact steady state 
probability distribution Pst oc e~*=* was found by Kwon 
et al. [11] with 

^st^lq'^ -U-q , (15) 

where the symmetric matrix U ~ is given by 

U = {D + Q)-^F (16) 

with the anti-symmetric matrix Q = —Q^ satisfying 

FQ + QF^ ^ FD- DF^ . (17) 

The equilibrium detailed balance condition, Eq. ([B]), 
yields Q = and U = D~^F, which satisfies Eq. (g]) 
as expected. Thus the existence of nonzero Q implies 
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the breaking of the detailed balance and results in the 
non-zero steady state current as jst = — (Q • ^^st)Pst- 
The general solution for Q can be found in a series form 
by using the Jordan transformation for asymmetric F or 
in an integral form in the frequency space [26| . 

The time-dependent PDF solution for the Fokker- 
Planck equation of Eq. ^ can be formally found from 
the path integral, using Eq. ([5]), as 

P{q,t;l)^ f dq{0)P{q{0)) f D[q]e~ ^^'^'^^^^^-'^ , 



(18) 

where P(g(0)) is the initial PDF and the Lagrangian L 
reads as 

L{q,q) ^\iq + F ■ qf ■ D-^ ■ {q + F ■ q) . (19) 

/ D[q] - ■ ■ denotes the integration over all paths reaching 
a fixed final state q{t) at time t, starting from an initial 
state 9(0), with proper normalizations. The source field 
1(t) is introduced for later use to generate time correla- 
tion functions and moments of q. 

We calculate explicitly the time-dependent PDF 
P{q,t) = P{q,t;l = 0) with the initial Gaussian PDF 
of P(g(0)) cx e-^'?'^(0)-'^('')-9(0) ^^^Yi a symmetric ^(0). 
The detailed calculation steps are given in Appendix [Al 
From Eq. (|A27p . we get 



P{q,t) = \dei{2TTA-\t))\-^/^e 



(20) 



where the symmetric matrix A ^{t) is given from 
Eq. (|X23l) as 



A-'{t) = +e-^\A-\0) - U-')e 



(21) 



As Vmvt^ao A~^{t) = U^^ for positive definite P, the 
steady-state solution of Eq. (fTSi) is recovered. 

The differential equation for A~^ can be derived from 
the recursion relation in Eq. (|A18p as [s^ 



dA- 



dt 



2D - FA- 



A~'F 



1 l?T 



(22) 



One can easily show that A~^{t) in Eq. ([21]) is the solu- 
tion of this differential equation. 

Now we define the generating functional for time- 
correlation functions and cumulants of q as 



Z[1{t)] = / dqP{q,t;l{T)) 



In Appendix |B1 we find from Eqs. (|B6|) and (|B7 



(23) 



(24) 



(r)e-(^'-")^ , r'>r, 



where r(r, r') = r^(T',T). Then we can compute the 
time average of any functional A of path ^(t): 



{A[q{T)]) = A 



' S ' 




[Sl{r)\ 


Z[l] 



(25) 



with Z[Q] = 1. For example, we get the two-time corre- 
lation function as 



('Za(r)g,(r')) 



5Ut) Slp{T') 



Z\l] 



1^0 



= T^piT,T')=Tp^{T',T) . (26) 

IV. NON-EQUILIBRIUM WORK PRODUCTION 

The work production by the NEQ force g along the 
path (/(r) is given by Eq. ([TU| with ^ = 0. As discussed 
in Sec. II, we have some arbitrariness in choosing the 
energy functional $((j), thus also the NEQ force q{q) in 
Eq. 0. 

In the case of the linear drift force with D^^ ■ f = 
-D~^F ■ q, Eq. ^ becomes 



(27) 



If D~^F is symmetric, the detailed balance condition, 
Eq. ^ is satisfied and one may choose $ = ^q^ -{D'^F)- 
q with g ~ 0. So we get no NEQ work production with 
this choice of the energy function, as expected. 

When D~^F is not symmetric, we must have a nonzero 
NEQ force g. The energy function $(g) can be written 
in general as 



Hq) ^^q^-Gs-q. 



(28) 



with a symmetric matrix Gs = which is a part of 
D~^F. Then we can divide D~^F into the symmetric 
part Gs and the remainder Ga- 

D-^F^Gs + Ga, (29) 

and the NEQ driving force is given as 

g{q) = -Ga-q. (30) 

There is no unique way to determine Gs or Ga out 
of D~^F. One natural possible choice is to enforce Ga 
anti-symmetric {Ga = —Ga), such as 

Gs = ^{D-^F + F^D-^) = Gs , 

Ga^^{D-^F-F^D-^) = Ga , (31) 

which will be called the anti- symmetric (AS) choice. 

Another interesting choice, called as the steady-state 
(SS) choice, is 



Gs = U , Ga = D~^QU 



(32) 



1^0 



If we take an initial Boltzmann distribution with this 
energy function, the system stays in the NESS from the 
beginning. 
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In general, one can choose 



where 



Gs — Gs + 5Gs , 

Ga = Ga — SGs , 



(33) 



with an arbitrary symmetric matrix 6Gs ■ From Eq. (|10p , 
the NEQ work production during time t is given by 



W[q] 



drf • Ga • g = W + (5W , (34) 



where W is the NEQ work production in the AS choice 



W[q] = - / drf -Ga-q 
Jo 



(35) 



Contribution 6W from the additive symmetric matrix 
SGs is given by 

6W= [ drf ■ SGs ■ q 
Jo 

= ^q'^it) ■ SGs ■ q{t) - ■ SGs ■ qiO) , (36) 

which comes only from boundaries and also exactly com- 
pensates the additional energy term due to SGs in the 
energy function $. 

It is important to note that the heat Q = A<1> — W is 
independent of the choice of SGs , in contrast to the NEQ 
work W. In the long-time limit, SW becomes negligible 
as the NEQ work usually increases incessantly in time. 
Thus the main contribution to the NEQ work production 
in the steady state comes from the purely anti-symmetric 
part VV. 

We now consider the generating function ^(A) for the 
PDF of the NEQ work production P(W) as 

GiX) = (e-^^) 

dq{t)dq{0)P{q{0)) f i^[q]e- /o ''-^(^^^-^WM ^ 



dWP{W)e 



-AW 



(37) 



with the initial equilibrium Boltzmann distribution 
P(g(0)) cx e-59^(o) Ga-9(o). The PDF of the work pro- 
duction P(yV) can be obtained formally by 



(38) 



We calculate G{X) in a similar way in which the path in- 
tegral is computed in Appendix |X] for the time-dependent 
PDF, P{q,t). As W is quadratic in q, the integral in 
Eq. ([37|) is basically the same as the integral in Eq. (fT8|) 
except for the final integral over q{t) with the modified 
Lagrangian as 

L ^ ^iq + F ■ qf ■ D-^ ■ {q + F ■ q) + ^q^ ■ A - q , (39) 



F = F- 2\DGa 
A = i [f^D^^F - F^D-'^F 

= X (F^Ga - GaF + 2X GaDGa) 



(40) 
(41) 



The contribution from only modifies the initial 

and final distribution according to Eq. ([36]). 

Before going further, we briefly comment on the 
discrete-time representation of the path integral and the 
work W. In Appendix]^ we perform the path integral in 
the discrete-time representation. Choice of the qir) value 
between the discrete time interval At does not affect the 
PDF at final time t in the limit of At — >■ 0. However, 
it is well known that one should choose the midpoint q 
value for the definition of the work for the correct de- 
scription [3l[. So VV in Eq. ([55)1 should be written as 



1 ^ 



i=l 



N 



-^{Qi - Qi-l)'^ ■ Ga ■ qi-1 



(42) 



1=1 



where = q{ti) tor i = 0,...,N with At = t/N. Note 
that it is first expressed in the mid-point representation 
but becomes identical to the so-called pre-point repre- 
sentation due to the anti-symmetricity of Ga- In Ap- 
pendix]^ all calculations are done in the pre-point rep- 
resentation for convenience. 

First, we perform the path integral of Eq. (|57|) with- 
out the final integral over q{t) . The integration procedure 
is basically identical to the case for the time-dependent 
PDF calculation in Appendix except for the differ- 
ent initial condition and the modified Lagrangian. Let 
A{t;X) be the modified kernel for A{t) of Eq. In 
the discrete-time representation, the recursion relation 
Eq. (jAlSp in Appendix \X\ is modified as 



2AtD + V (A^ 



AtAj 



(43) 



with V = 1-AtF. Taking At ^ limit, the differential 
equation for A~^ can be derived, 



dA- 



dt 



= 2D- FA-^ - A-^F'^ - A~^AA~^ 



(44) 



In contrast to Eq. (P^ . this is a nonlinear differential 
matrix equation, which can not be solved analytically 
in general. Using dA^^/dt = —A^^{dA/dt)A~^, we can 
rewrite this equation as 



dA 

dt 



-2ADA + AF + F' A + A 



(45) 



The initial condition is given as ^(0; A) = G'^ + {l—X)SGs, 
where the A-dependent term comes from SW. We solve 
this equation numerically for a specific case in Sec. VI. 
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We obtain 5(A) after integrating over the final q{t) 
(Jtv in Eq. ([571). leading to 



SiX) 



det(G. + 6Gs) 



-1/2 



N-1 

n 



det(A, + At A) 



det{A 



i+l/ 



-1/2 



(46) 

Using Ai^i ~ Ai + {At)dAi/dt for the denominator, we 
finally get 



Ing(A) 



1 



drTr A 



di(T; A) 



A-\t;X) 



(47) 



^2 " i det(G, +(5G,) 



drTiiAD - Z') 



_ 1 / det(i(^;A) + A(5G,) \ 
'2""^ det{G,+SGs) J 

It is not possible to perform the integral of Eg. ([55]) to 
find P{W) in a closed form. However, the explicit form 
of G{X) reveals many interesting properties of P{W). For 
example, as \nQ is not quadratic in A, P{W) is not Gaus- 
sian in general. Due to the logarithmic boundary term, 
the divergence Q{X) may appear in A, which determines 
the asymptotic behavior of a non-Gaussian tail of P(>V) 
for large |yV|. This will be investigated more in detail in 
Sec. VI. 

The fluctuation theorem yields 



g{x) = g{i - A) 



(48) 



by substituting A with e~^^ in Eq. ([13]). The work 
yV[q] in the reverse path should be the same as — ^^[(j] in 
the forward path, and the forward and reverse path are 
identical with the same initial conditions with the same 
energy function. It seems not easy to prove Eq. pS)) for 
general A, directly from Eq. (|T7)) . 

However, we can prove tj(l) = (e~^) = 1 easily, 
which corresponds to the Jarzynski equality for a time- 
dependent potential. For A = 1, A{0) ~ Gs, F = 
DGs - DGa, and A = GsDGa - GaDGs- We can show 
the initial state is the fixed point of Eq. p4)) or (|45|). 
i.e., A{t) = Gs. Then the logarithmic part vanishes in 



Eq. (gT]). We can also see Tt{AD 
Hence g{l) = 1. 



F) = TiDGa = 0. 



V. CUMULANTS OF NEQ WORK 
PRODUCTION 



In this section, we calculate the cumulants of the work 
production by using the two-time correlation function 
r(r,r') in Eqs. ([24]) and ((Ml). For simplicity, we take 
the AS_ choice where G^ = Ga , SGs = 0, SW = 0, and 
W = W. However, in the long-time limit, all results are 
choice- independent. 



First, consider the first cumulant of W in the discrete- 
time representation, Eq. (|^^ as 



N 



i=l 

N 

= ^Trr,,,_iGa 

i=l 
N 

= Y^Tiil - AtF)A-_\Ga 

i=l 

N 

= -Y,{At)TTFA-_\Ga 



i=l 

ft 



- / dTTTA-\T)GaF 

Jo 



(49) 



where we use qJ^i-Ga-qi-i = 0, T^.i^i = {I- AtF)A^\, 
and Eq. (gH). 

In the long time limit, we can replace A~^ by U^^ and 
find 



(W) ^ -t TiU-^GaF 
= -t TlQGa 
= t TtQF'^D'^ , 



(50) 



which should be choice-independent. Since the energy 
difference (A$) is finite, (W) ~ — (Q) measures the en- 
tropy production piled up in the reservoir with the mean 
rate 



(a) = t-\W) = TtQF^D~^ 



(51) 



It is expected to be positive from the fiuctuation the- 
orem, which implies that the second law of thermody- 
namics should hold for general NEQ phenomena with no 
thermodynamic origin. 

The second cumulant of W can also be found as 

(W^), = (W^) - (>V)2 

= Y,{{qf -Ga-m-i) {qf ■ Ga ■ q,^i)) 
I 

= Tr {—TijGaTj-ii-iGa + Tij-iGaT ji^iGa) (52) 
where the summation Y^' • = 2 y^,;^,, + y^,_, due to the 
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symmetry between i and j. Using Eq. (I24[) . we can write 

= E 2^Tr{-e-^(*'"*^)ATiG„^7_\e-^'"(*.-i-*.-i)G, 

i I 3 = 1 



-Tr [-AT'GaAT^Ga 



In At limit, we get 

-2Tr /" dTe-''^''{F^Ga-GaF)Fe 
Jo 



(53) 



From Eq. ([55]) . (cr^)c ^ It qualitatively agrees 

with the fluctuation theorem for the entropy production. 
However, the fluctuation theorem predicts s = 1, while it 
seems not equal to 1 if estimated by assuming the Gaus- 
sian distribution. This implies that the PDF for the work 
and entropy production is in general non-Gaussian and 
non-Gaussian tails make a significant contribution to the 
exact theorem. 

More information on P{yV) might come from higher 
cumulants in W. In principle, it can be done system- 
atically by using Eq. ([25|) . but it is too complicated to 
proceed further calculation in detail. 



-Ft 



X / dr'e^^' A-\T')GaA-\T')e^^^' 



-2Tr / dre-^ ^{F^Ga-GaF)e 

X r dT'e^^'A~\T')GaA-\T')e^"^' 



-2Tr / dTGaFA-'{T)GaA-'iT) 
Jo 

-^Tr {GaA-\t)GaA-\t) - GaA-^{0)GaA-'{0)){^4) 



VI. EXAMPLE: DIFFUSION IN TWO 
DIMENSIONS 

Now, we take an example of a two dimensional diffu- 
sive motion, q = {x, y)^ for more explicit calculations. 
Consider 



(59) 



By the orthogonal coordinate transformation, D can be 
diagonalized and the calculation goes simpler. However, 
we keep the present form of D in order to examine the 
effect of noise correlations. If FD — DF^ ^ 0, we have 
a nonzero anti-symmetric matrix Q, a measure for NEQ, 
which can be obtained easily from Eq. ([TT)) . 



Ft 




( fci 














fc2 ) 







where A~'^(t) is given by Eq. 

Using the identity of Eq. (|A21[) found in the Ap- 
pendix one can perform the above integral in prin- 
ciple. In this paper, rather than reporting the exact 
time-dependence of Eq. ([M)) . we compute the long-time 
behavior by keeping only the most dominant contribu- 
tions, 

^ -t Tr (C + E) {F^Ga - GaF) 

-2t TlGaFU^^GaU^^ 

= 2t Tr [GaF{E - G)] (55) 

where the matrix G is anti-symmetric, defined by 

C = U-^GaU-^ , (56) 

and the matrix E is symmetric, determined by 

EE + EE'^ = EG - GF'^ . (57) 

Note that {{t-^W)^)c ^ i"^ It implies that the PDF 
of the entropy production rate, a = t^^W, shows a sharp 
distribution with the mean value (cr) found in Eq. (j5ip 
and the variance of order t~^. Assuming it as Gaussian, 

P{cr) ^ e 2<''^>c , we obtain 



q 

-q 



q 



e(fci - fc2) -I- 7Ki — aK2 
ki + k2 



(60) 



The system goes to equilibrium for ^ = 0. The conven- 
tional Gibbs-Boltzmann (GB) type equilibrium is a triv- 
ial case where the force is conservative (ki = K2), and the 
noises are identical and independent (a = 7, e = 0). The 
equilibrium PDF is given as e^" £ix,y)^ ^^^^^ £{x,y) = 
^{kix^ +2Kixy+k2y^) ■ There are also non-trivial equilib- 
ria possible even for a non-conservative force (ki 7^ K2) 
and non-identical/correlated noises (a 7^ 7, e 7^ 0) as 
long as g = 0. In this case, the equilibrium PDF is given 
as e~^' ^■'J. There seems no fundamental difference 
among many equilibria from the view point of our pa- 
per in the sense that they are all preserving the detailed 
balance. 

Now we consider the NESS for q ^ 0. From Eq. ([CT]) . 
the entropy production rate in the long time limit is given 
as 



TtF 
detTo' 



(61) 



where dct D = a-f — and TrF = ki + k2- 

The second cumulant of the work can be obtained from 
Eq. (|55|) . The matrices G and E are found from Eqs. ([55)) 
and (|571) as 



^(-) _^ts. ^ 2(a) 



Pi-a) 



(58) 



G = 



c 

-c 



E = 



^2ki ki — ^2 
TrF V fci - ^2 2k2 



, (62) 
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'^ith ^' = 2dctJdctf g- Then, (cr^) is found as 

For the stabihty of the NESS, we assumed at the be- 
ginning that the matrices D and F are positive definite. 
Thus, dcti? > and TrF > 0, which guarantee the 
positivity of (yV)c and (W^)c- From Eq. ([SHI), it can 
be shown that s — dot Z?/ (dot Z) + q^) < 1, which indi- 
cates that P(cr) is more distributed than the Gaussian 
distribution for nonzero q. We observe that aU higher- 
order cumulants are also of order t, i.e. (yV")c ^ t, so 

We have performed numerical analysis to confirm our 
analytic results and gain more insights on the work distri- 
bution function P(yV). Here we present numerical data 
taken at (^1,^2) = (4,1), (^1,^2) = (2,1), (0,7) = (1,1), 
and e = sinO with 6 = 0.1. This system has a NESS with 
a nonzero value of q. We expect that the results do not 
depend on a specific choice of parameter values as far as 
q is nonzero. 

For Gs and Ga, we adopt the AS choice given in 
Eq. (|5T|) for convenience. Hence, the system is assumed 
to have an initial probability distribution 

P{q{0)) = |dct(27rA-i(0))ri/2g-ig-(o).A(o).,(0) 

with A{0) = Gs and the work for a path g(r) is obtained 
using Eq. 

The generating function t?(A) defined in Eq. ((37)) 
can be estimated by direct numerical integrations of 
the Langevin equation in Eq. ([1]). One starts from 
an initial state q(0) drawn from the initial distribution 
P{q(0)). The Langevin equation is then integrated with 
discrctizcd time intervals 

q{T + At) = q{T) - F ■ q{T)M + AW{t) 

where AW^(t) = {/S.Wx{t)t /S.Wy{T))'^ arc correlated ran- 
dom variables constructed as 

\ AWy{T) j-^^^^y sinf? cose )\ rjyir) J ' 

Here, rj^^y^r) are independent and identically distributed 
Gaussian random variables with zero mean and unit 
variance. One can check easily that such random vari- 
ables AW{t) satisfy the required correlation property 
{AW{t)AW{t')'^) = 2{At)DS{T - t'). 
The work production is estimated as 

W(t + At) = W(t) - {q{T + At) - q{T))^ ■ Ga ■ q{T) . 

Repeating Ns independent simulations, one obtains a nu- 
merical estimate 

5Nit;X) = ^f:e-'^^(^K (64) 



We can also utilize our analytic expression for Q in 
Eq. (|T7)) in order to get a more precise numerical esti- 
mate. We first solve the NLDE for A{t; A) in Eq. gS]) 
with the initial condition A{0) = Gs- The solution is ob- 
tained numerically in discretizcd times using the recur- 
sion relation A(r+At) — A{T)+At '^^[^^ . The generating 
function Q can be evaluated easily using the numerical so- 
lution A. The generating function evaluated numerically 
using the analytic expression will be denoted as Qa- 
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FIG. 1. Comparison of C/]v(A) and Ga{^) at i = 0.5 in (a) and 
t = 2.0 in (b). Symbols represent Qm and lines represent Qa- 
Location of the symbol A represents the Jarzynski equality 
5(1) = 1. In (a), Qa diverges continuously as A approaches a 
threshold. On the other hand, in (b) , Qa remains finite up to 
a threshold and diverges discontinuously beyond it. 

In Fig.[Tl we compare the results GN{t; A) and GA{t', A) 
at t = 0.5 and 2.0 obtained with At = 0.0001 and 
< 10^- The two methods yield almost identical 
results for small values of |A|. Both data confirm the 
Jarzynski equality Q{X = 1) = 1 and the Crooks fluctua- 
tion theorem C/(A) = Q{l — X). However, there is a notice- 
able discrepancy at larger values of |A|. Even the Crooks 
fluctuation theorem seems to be violated in the data. 
One might suspect a finite At as a source of systematic 
errors. We have also taken data with A< = 0.01,0.001, 
and 0.0001 and found no significant difference, which 
means that At = 0.0001 is already small enough. In 
fact, the discrepancy is due to limited sampling in ob- 
taining Qn- When |A| is large, Q{X) is dominated by rare 
events with large |yV|. If one compares Gn obtained from 
Ns = 10^, 10^, and 10'' samples, there are strong fluctu- 
ations at large values of |A|. This means that the tail 
property of the work distribution function P(yV) cannot 
be accessed from numerical simulations even with 10^ 
samples. On the contrary, the analytic formalism allows 
us to study the work distribution function in detail with- 
out any statistics problem. 

Figure [1] shows that G{t; A) becomes singular at A = Aq 
and 1 — Ao with a t-dependent threshold Aq = Xo{t) > 
1. The singular behavior is evident in Fig. [2ja), where 
we plot \nG{t]X) as a function of t at several values of 
A > 1/2. It suffices to consider A > 1/2 because of the 
symmetry ^;(A) = ^^(1 - A). When A < A^ ~ 1.962(1), 



FIG. 2. (a) Plot of lnQ{t; A) against t at several values of A. 
(b) Threshold curve Xo{t) above which Q{t; A) is infinite. The 
symbol represents a point {tc,\c) — (0.86, 1.962). 



Q(t;X) remains finite for all t. On the other hand, it 
diverges at t — tc ~ 0.86(1) at A = Ac, and diverges 
a.t t < tc when A > Ac. From these plots, we conclude 
that G{t; A), being viewed as a function of A, diverges at 
i-dependent thresholds A = Xa{t) and A = 1 — Ao(i). The 
threshold, numerically determined, is drawn in Fig.[2][b). 
Figure 2 also allows us to conclude that Q{t; A) diverges 
continuously as A — >■ Xo{t) (see the solid line in Fig.[2jb)) 
when t < tc- On the other hand, when t > tc, G{t;X) 
displays a discontinuous jump to infinity at A = Ac (see 
the dashed line in Fig. [UJb)). Interestingly, Fig. [2jb) 
resembles a phase diagram of a system having a tricritical 
point where a continuous phase transition line turns into 
a discontinuous phase transition line. 

Origin and nature of the divergence are understood 
from the analytic expression for \nQ in Eq. (jlT]) . Due 
to the logarithmic boundary term, G(t] A) is well-defined 
only when det{A{t'; X) + XSGs) is positive for all t' < t. In 
contrast, there is no singularity in the bulk term for any 
t. From the numerical solution of the NLDE, Eq. (|45|) . we 
observed that dct{A{t; A)) {6Gs ~ with the AS choice) 
behaves as 

det(i(t; A)) ~ a(Ac - A) + bit - tcf (65) 

near A = Ac and t — tc with positive constants a and 6, 
see Fig. [31 This behavior explains the singularity in Q. 

When t < tc, the determinant becomes zero at A = Ao 
with 

Ao Ac + -(t-^c)^ . (66) 
a 

Consequently, Q diverges continuously as 

g(t;A)~(Ao-A)-i/2 (67) 

as A approaches Aq from below. 

The determinant is always positive when A < Ac, while 
it becomes negative at t < tc when A > Ac. Hence, when 
t > tc G{t;X) remains finite up to A < Ac and diverges 



FIG. 3. Time evolution of det A{t; A) near t = tc and A — Ac. 
When A < Ac — 1.962, the determinant seems to be always 
positive for any t. In contrast, it becomes negative at finite 
t for A > Ac and the time evolution afterward (the dashed 
line) is meaningless. When A = Ac, it is tangential to the x 
axis (the horizontal dotted line) at t = tc (the vertical dotted 
line). 

discontinuously at A = Ac. As A approaches Ac from 
below, it behaves regularly as 

g{t;X) ^ h + a{Xc- Xy (68) 

with the t-dependent constant h. 

The singularities in ^^(A) = / dWP{W)e-^^ at A = 
Ao and 1 — Aq indicate that the PDF P(yV) has exponen- 
tial tails 

r W-'' e-^/^+ for W ^ oo 
P{W) - <^ (69) 
I i-Wy e-^l^- for W ^ -oo 

with characteristic works yV-|_ (> 0) and yV_ (< 0), and 
possible power-law corrections with exponent r. The 
power-law prefactor is necessary in order to account for 
the way how Q becomes singular at A = Aq. From the 
Crooks fluctuation theorem, P(W) = e^P(-W), the ex- 
ponent r should be the same for both tails, and it suffices 
to consider one of the tails. 

It is easy to check that the negative tail yields ^(A) ^ 
(1/|W_| - Xf-^. Comparing it with Eqs. (gZl) and 
we find that the characteristic work is given by W_ = 
— I/Aq and that the exponent is given by 

r i for t < tc , 

(70) 

\2 ior t>tc . 

The positive tail has yV+ = l/(Ao — 1) and the same 
exponent r. 

Numerical data are consistent with the tail property in 
Eq. The PDF P{W) was obtained from numerical 

simulations of Ng = 10'' samples. Figure Hl^a) presents 
the plot of P(yV) against W at several values of t in 
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FIG. 4. (Color online) (a) Plot of P(W) versus W at several 
values of t. (b) Plot of P(>V)e'^o-i)^ versus W for W > 0. 
The solid and dashed lines are guides to eyes with slope —1/2 
and —2, respectively. 



the semi-log scale. As expected, P(W) becomes more 
distributed and the mean work production (W) increases 
with time t. Moreover, the exponential tails are clearly 
seen in all plots with characteristic works (slopes of plots 
in Fig. m^a)) saturating with time t. 

We also present the log-log plot of P{W)e^^°^^'^^ 
against W for > in Fig. gljb). We use the 
threshold value Ao(i) obtained from the singularities in 
Qai which are shown in Fig. [2jb). One can see that 
the P(yV)e'^^°~^^^ has a power-law tail as predicted in 
Eq. dini)- When t = 0.2 and 0.5 (< tc), the power- 
law tail is manifest and the exponent is in good agree- 
ment with the analytic prediction r = 1/2. When 
t = 2.0 and 5.0 (> tc), the power-law sets in at larger 
values of W and the exponent value is drifting with in- 
creasing W. In this case, we need much more samples 
[Ns S> 10^) to get good statistics for rare events at large 
W, in order to extract reliable quantitative information 
on the power-law tail. Nevertheless, one can see that 
the exponent value becomes close to ?■ = 2 for large t. 
Pd t = 1.0, the crossover effect dominates the numerical 
data, since t is too close to tc — 0.86. 

It is not surprising to see the exponential tail with a 
power-law prefactor with the exponent r = 1/2. This ex- 
ponential tail has been also observed in other NEQ sys- 
tems [H-dil. When the PDF is given as P{q) - e"**'' 
with the energy-like function $((7) quadratic in q as in 
our case, the PDF for any quantity A also quadratic in 
q should become P{A) ^ A^^l'^e^^l-^" with the charac- 
teristic value of Ac- As the work W is simply quadratic 
in q and proportional to t, at least for a short time i, 
i.e. W ~ tq^ as in our case, our finding in Eqs. (|69|) and 
([70)1 can be understood with Ac increasing in time t. 

However, there is a sharp dynamic phase transition 
aX, t = tc beyond which Ac is constant of time and 
the power-law exponent changes from 1/2 to 2. This 
implies that the characteristic positive work production 
W+ = l/(Ao(t) — 1) increases monotonically in time for 
t < tc, but saturates to a finite value l/(Ac— 1) at t = tc- 



For t > tc-, the characteristic work production remains 
unchanged and only the power-law prefactor adjusts the 
PDF accordingly. We do not have an intuitive under- 
standing for the dynamic phase transition at this mo- 
ment, except that the transition time tc may be related 
to the intrinsic relaxation time of the system. It calls 
for a further study to understand this phase transition, 
which is currently under investigation. 



VII. SUMMARY AND DISCUSSION 

We have studied NEQ fluctuations for high- 
dimensional diffusion dynamics driven by a linear drift 
force. The drift force is not derivable from a scalar po- 
tential function and the noises are not identical white 
noises in components with possible correlations. In gen- 
eral, these NEQ features generate a circulating probabil- 
ity current at the steady state, with the drift force inward 
to the origin (the force matrix F is positive definite) . It is 
interesting to notice that the equilibrium can be restored 
with a specific combination of these two NEQ features by 
satisfying the detailed balance condition {DF = F^D). 

Recent experiments on an optical trap report that par- 
ticles can be confined in a field-induced potential well, 
approximated by an asymmetric harmonic potential. It 
may be possible to apply our study to such experiments. 
However, there are some technical issues to be resolved 
in experiments. In particular, in order to break the de- 
tailed balance, the noises should be applied in such a way 
that their principal axes do not coincide with those of the 
potential well [25| . 

In this paper, we started with the Langevin equation 
with an arbitrary drift force and arbitrary additive noises 
in high dimensions without any thermodynamic origin. 
The generalized thermodynamic quantities like energy, 
work, and heat are defined, with which the Crooks fluc- 
tuation theorems are reproduced. For the case of the 
linear drift force, we derive the exact evolution func- 
tion of the PDF analytically. More importantly, we ana- 
lyzed the work PDF, P{W) through the generating func- 
tion method, which showed an interesting dynamic phase 
transition in the exponential tail shape of P{W). As the 
tail property is governed by rare events, it is crucial to 
develop an analytic theory to show this subtle dynamic 
phase transition, which can be hardly identifled by nu- 
merical simulations only. Even though we expect that 
this phase transition is one of the generic features found 
in the NEQ systems described by Langevin equations, it 
would be a big challenge to understand analytically how 
and when this can arise. 

Finally, in the case of nonlinear forces, it has been 
found that there is an additional current which shifts the 
probability maximum from the fixed point at the origin 
[26|. It may have the same origin as the noise-driven 
directed current in extended systems, [13, [l^. The in- 
gredients for those currents are: (a) non-linearity of the 
inward drift force and high-dimensional noises for the 
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confined system (b) asymmetric (ratchet-type) potential 
and an additional noise for the extended system. The 
NEQ fluctuation theory for this non-linear diffusion dy- 
namics will be another challenging topic. 
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Appendix A: Time-dependent PDF 

We keep the source field for later use to compute the 
generating functional. The initial probability distribu- 
tion is chosen as 



P(<Z(0)) =AAoe-39^(o)-A(o).9(o) ^ 



(Al) 



with AAq = |dct(27rA-i(0))|-i/2 ^nd A^{0) = A{0). We 
consider the sequence of discrete times, = to < < 
t2 < ■ ■ ■ < In = t with time interval At = t/N in the 
N ^ oo limit. Denoting qi = q{ti) with qq = q{0) and 
Qn = q{t) = 9, and Aq = A{0), we write P{q,t;l) in the 
pre-point representation (set (/(r) = qi_i for ti-i < r < 
ti) as 



» N-l 

P{q,t;l)^No J n 



dqi 



^ |det(47rAti:')|i/2 



'7E£i(9.-9.-l+AtF-q,_i)^-D-i-(g,-gi_i+Ati^-gi_i) 



xe 



» N-l 

A/'o/n 



dqi 



|det(47rAii:))|i/2 



N 



(A2) 



The transfer matrix r,;.i_i is given as 



lnri,i-i(/j-i)= {qi-Vq,^i)'^ ■ D ^ ■ {qi-Vq,^i) 



Integrating over we obtain 

° J |det(4^Aii?)|i/2^i^o'' 
= |dct(i)Bo)r^/^e-^9r-^i-'?i 



where 



1 

2At 

Bo = V'^D-'^V + 2AtAo , 



Ai = Tr^ (D-^ - D-^VB^^V^D-^) (A6) 



(A7) 



where Ai and Bq are both symmetric. Then we can find 



T2,i(Zi)/o = T2,i([i)e-^«"-^-^^ 
x|det(7^i?o)r'/'e('^*'''"-^o-^-'o 



(A8) 



where li = h+D ^VB^^ ^-lo- Notice that the integration 
of Ti.o over qo results in a change h in the source field of 
T2_i. It happens iterativcly for subsequent integrations 
over the rest of qi. Wc show it explicitly for the next 
step, 

^ J |dct(47rAtI?)|i/2^2,il'iJe 



|det(i:»Bi)ri/2e-|9J'-^2-92 
xe 



{Atfl^-B^^-li + {At)l'^-B^W^D-^-qi 



(A9) 



where 



1 



A2^Trx7 (^"' - D-'VB^W^D"') (AlO) 



2At 

Bi = V'^D-^V + 2AtAi . 
Then we get 

^3,2(^2)/! =T3,2(^2)e'^ 



(All) 



Ql ■A2-qi 



x\det{DBi) 



-1/2 {Atfl^-B-^-h 



(A12) 



where I2 = h + D ^VB^ ^ -h- We repeat the integrations 
over all qi up to i = — 1. Finally, we reach the result: 



N-l 



where 



(At)Zf_i -q,.!, 



V = 1- AtF. 



(A3) 
(A4) 



PiqN,tN;l) = |det(27rA^i)|-i/2 Yl |det(i^B,)|-i/2 

1=0 

We have the recursion relations: 

= ^^{d-' D-WB-_\V^D-') (A14) 

= V'^D-W + 2AtA,_i (A15) 

~k = k + D-^VB,^i-~k_i, (A16) 

for i = 1, . . . , A^ with Iq — Iq. 
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Ai determines the intermediate probability distribu- 
tion function at time ti. We can find a simple recursion 
relation for A~^. First we observe the following relation: 

A, = ^ f^D-' -(^D + 2AtD{V^y^A,^,V-'Dy' 

= {V'^)-^A,_iV-^d(^D + 2AtD{V'^y^A,^iV-^Dj 
= (^2AtI? + l/A7_\F^) \ (A17) 
Therefore we get 

A-^ = 2AtD + VA-\V'^ , (A18) 

which leads to 

A-^ = 2AtD + 2AtVDV^ + V'^A~}^{V'^f 

^2At^VW{V^y +V'A-^^{V'^y . (A19) 
In the continuum limit, we have 

A-\t)=2 f dre-^-De-^"- + e-^'A~\Q)e-^"' . 
Jo 

(A20) 

We use the integral identity as 



which simplifies the normalization factor in Eq. (jA13p as 



N~l 



/2 



I det(27rAo 1)1-1/2 \[ |det(7^i?,)|-i 



A det^o det Ai detAjv-i 

det(27r^Q ) 



det Ai det A2 det A 



N 



-1/2 



|det(2^A^i)ri/2 



(A26) 



dTe-"''Ce 



-H- 



(C + E) - e-^'^C + E)e'" * , (A21) 



for arbitrary matrices C and H where E is determined 
from the equation 



HE + EH^ ^ HC - CH^ 



(A22) 



Note that E is antisymmetric (symmetric) if C is sym- 
metric (antisymmetric) . 

Here, with C = D and H = F, the antisymmet- 
ric matrix E becomes equal to Q from Eq. ([TT]). Then 
H-^{C+E) = F-\D + Q) = C/-1 fromEq. There- 
fore, we have 



A-\t) = [/-I + e-^*{A-^{Q) - [/-i)e- 



t 



(A23) 



Using Eq. (|A18|) . the recursion relation for S^-i, 
Eq. (|Al5p . becomes 

= Ai_iy-i {yA-\v'^ + 2AtD^D-'^V 
= Ai_iV-^A-^D-W . (A24) 
Therefore we get 

det Ai_i 



det(I?B,_i) 



It is the correct normalization factor for the final time- 
dependent probability distribution, P((7, t) = P{q, t; I = 
0): 

P(g,t) = |det(27rA-i(i))ri/2e-^«^-^W-« . (A27) 



Appendix B: Generating functional 

The generating functional Z[l] = J dqP{q,t;l) is ob- 
tained by integrating out Eq. (|A13p over qpf-. 



1 

= £2 



(Bl) 



where = 2{At)Aiy. Investigating the recursion rela- 
tion for k, Eq. (|JT6| . we first show that for i > j 

Tu = 2AtSri + r,+i,,+iQ (B2) 



Til — TiiCi—iCi-2 ■ ■ ■ Cj 



(B3) 



where 



a - D-^VBy^ - A,+^VA-^ . (B4) 
Using Tnn ^ ^Tv^ and Eq. (|A24p . 

rN-i,N-i = 2AtV-'DANVA]^y, + C^_iA^iCjv-i 

= 2AtV-'DANVA]^l^ + A^i_,F^AjvFA^i_i 

(-1 T/T\ a ta 4-1 



= V-'[2AtD + VAyy_^V' jA^VAy^l^ 
= A 



We can prove by induction 



By noting 



we find 



■ ■ ■ c, = A^v'-' A-^ 



= V'-'AJ^ 



(B5) 
(B6) 



(B7) 



Note that for j > i, Tij ~ F^. In the continuum limit. 



det A, 



(A25) wc finally get Eq. (HH) 
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